Non-equilibrium clustering of self-propelled rods 
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Motivated by aggregation phenomena in gliding bacteria, we study collective motion in a two- 
dimensional model of active, self-propelled rods interacting through volume exclusion. In simulations 
with individual particles, we find that particle clustering is facilitated by a sufficiently large packing 
fraction rj or length-to-width ratio k. The transition to clustering in simulations is well captured by 
a mean-field model for the cluster size distribution, which predicts that the transition values Hc of 
the aspect ratio for a fixed packing fraction r/ is given by Kc = C/rj — 1 where C is a constant. 

PACS numbers: 87.18.Bb, 05.70.Ln, 87.18.Ed, 87.18. Hf 



Introduction. — Emergent large-scale patterns of inter- 
acting self-driven motile elements are observed in a wide 
range of biological systems of different complexity: from 
human crowds, herds, bird flocks, and flsh schools Jj to 
multicellular aggregates, e.g. of bacteria and amoebae 
|2| as well as sperms |2J. A recurrent question is how 
these entities coordinate their behavior to form groups 
which move collectively. At a theoretical level, several 
qualitative approaches have been made to incorporate 
the diverse collective behaviors of such different systems 
in a common framework [ij, U, la] ■ More specific models 
for bacteria like E. coli as well as for amoebae like D. dis- 
coideum ^j, have been based on chemotaxis, a long-range 
cell interaction mechanism according to which individual 
cells move in response to chemical signals produced by 
all other cells. However, in some bacteria there is no 
evidence for chemotactic cues and cells coordinate their 
movement by cell-to-cell signalling mechanisms in which 
physical contact between bacteria is needed |g|. Con- 
sequently, one may ask how such bacteria aggregate in 
order to communicate. 

Another relevant aspect is the influence of the shape of 
the bacteria. The shape has been shown to be essential 
for individual motion of swimming bacteria 7] . In con- 
trast, the role of the cell shape for collective motion has 
remained mostly unexplored. It has been demonstrated 
experimentally [g] that migrating elongated amoeboid 
cells exhibit alignment effects similar to those reported 
in liquid crystals .Oj. A prominent example for collec- 
tive behavior with no apparent long range interactions 
are the striking patterns observed during the life-cycle 
of ghding myxobacteria, see e.g. [y, |l3- Earlier mod- 
eling work has reproduced many of these patterns in 
three dimensions assuming either perfect alignment [ll| 
or a phcnomcnological alignment force |l2l| . These mod- 
els have all considered patterns resulting from exchange 
of chemical signals, that are absent in an early stage of 
the myxobacterial life cycle. Nevertheless, a trend from 
initial independent motion towards formation of larger 
clusters of aligned bacteria is often observed ( Fig. 1). 



FIG. 1: Example for clustering of myxobacteria (M. xanthus) 
in the early stage of the life cycle, (a) Immediately after 
maturation of spores, (b) Afterwards, during the vegetative 
phase. Snapshots are taken from a movie (Ref. [13 (b)), the 
frame size is 40 x 30 fj,m^. Similar phenomena were seen in 
other bacterial species (c/. Ref. 0(a)). 



Here, we study a model of self-propelled rods that 
have only repulsive excluded volume interactions in two 
dimensions. We find that the interplay of rod geome- 
try, self-propulsion and repulsive short-range interaction 
is sufficient to facilitate aggregation into clusters. In 
simulations of an individual based model (IBM), clus- 
tering of self propelled particles (SPP) is observed for 
large enough packing fraction rj resp. aspect ratio k of 
the rods (see Fig. EJ. We define the onset of clustering 
by the transition from a unimodal to a bimodal clus- 
ter size distribution. A mean-field approximation (MFA) 
for the cluster size distribution is derived and reproduces 
the change from a unimodal to a bimodal shape upon 
increase of either rj oi k. The MFA yields a simple equa- 
tion Kc — C/?] — 1 for the critical rod aspect ratio, Kc, at 
the onset of clustering in line with the IBM simulation 
results, fitted with C = 1.46. If diffusion is added to 
the active motion (active Brownian rods), the clustering 
transition is shifted to higher values of k, whereas cluster- 
ing is absent for pure diffusive motion (Brownian rods) as 
well as for isotropic particles with k = 1. Hence, cluster- 
ing of particles with excluded volume interaction requires 
both active motion, i. e. a non-equilibrium system, and 
elongated particles (= rods). 

Individual-based model (IBM). — We consider N rod- 
like particles moving on a plane. Each particle is 







FIG. 2; Simulation snapshots of the steady states for different 
particle anisotropy k and the same packing fraction rj (a-c), 
and the same k and different 77 (d-f). Fixing -q — 0.24: (a) 
before the transition, k = 1; (b) almost at the transition, 
K = 5; (c) after the transition, k — 8. Fixing k = 6: (d) 
before the transition, rj = 0.18; (e) just crossing the transition, 
77 = 0.24; (f) after the transition, rj = 0.34. In all cases, 
particles A'^ — 100 and particle area a = 0.2. The arrows 
indicate the direction of motion of some of the clusters. 



equipped with a self-propelling force acting along the 
long axis of the particle. We assume that particles are 
submerged in a viscous medium. Velocity and angular 
velocity are proportional to the force and torque, corre- 
spondingly. The rod-shape of the particles requires three 
different friction coefficients which correspond to the re- 
sistance exerted by the medium when particles either ro- 
tate or move along their long and short axes. Inertial 
terms are neglected (overdamped motion) . Consequently 
the movement of the i-th rod is governed by the follow- 
ing equations for the velocity of its center of mass and 
angular velocity: 
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where w, , Wj^ refer to the velocities along the long and 
short axis of the rods, respectively, Q indicates the corre- 
sponding friction coefficients {Ce is related to the friction 
torque), C/^*-* refers to the energy of the interaction of 
the i-th rod with all other rods, and F is the magnitude 
of the self-propelling force. The motion of the center of 
mass x^*' = {vx , Vy ) of the i-th rod is given by 



;«=i;^'^cos0«-ff^,'^sin0« 



vf =vf^iAne^'^ -v{^cosi 



(2) 



Particles interact by ,,soft volume exclusion, i. e. by a 
potential that penalizes particle overlaps in the following 



way: 

t/W(xW,6l(*),x(^'),6l(^')) 

N 

= <^ E ((7-ao(x«,e«,x(^),0(j)))-^-7-'^) (3) 

where ao(x'^*\ 6''^*\x'^-'\ 6''^-''') is the area overlap of the 
rods i and j, 7 is a parameter which can be associated 
to the maximum compressibility, (3 controls the stiffness 
of the particle, and </> is the interaction strength. The 
simulations were performed placing N identical particles 
initially at random inside a box of area A with periodic 
boundary conditions. The values of the parameters are 
given in |l3| . 

There are three key parameters which control the dy- 
namics: i) persistence of particle motion, regulated by 
F, ii) the packing fraction 77, i.e., the area occupied by 
rods divided by the total area (77 = Na/A, where N is 
the number of particles in the system, a is the area of a 
single particle, and A is the total area of the box), and 
iii) the length-to-width aspect ratio k {k = L/W, where 
L is the length and W is the width of the rods). Simula- 
tions yield an increase of cluster formation with increas- 
ing K or 77, see Fig. El Individual clusters are defined 
by connected particles that have non-zero overlap area. 
Simulations can be characterized by the mean maximum 
cluster size, M, and the weighted cluster size distribu- 
tion, p{m), which indicates the probability of finding a 
given particle inside a cluster of mass 777,. Fig. 13^- shows 
that for a given 77, M seems to saturate after the critical 
Kc which is defined as the value of k for which the shape 
of p{m) changes from unimodal to bimodal. In Fig. ^p 
typical shapes of p{m) are shown: before clustering and 
corresponding to low values of k (circles), and after clus- 
tering and corresponding to large values of k (crosses). 
We define the onset of clustering by the emergence of a 
second peak in p{m). We have also tested the robustness 
of the model against fluctuations by inserting additive 
noise terms Rt/Ci in Eqs. Q, which correspond to a 
switch from active to active Brownian particles [a|. We 
found that clustering is still present in rods of the latter 
kind, albeit the transition is moved to larger values of 
K and 77. Clustering was absent in all simulations with 
purely Brownian rods {F — 0). 

Mean field approximation (MFA). — We have studied 
the clustering effects described above through a MFA by 
deriving kinetic equations for the number Uj of clusters 
of a given size j. The equations for rij contain terms for 
cluster fusion and fission. For the fusion terms we have 
adopted kinetic equations originally derived for coagula- 
tion of colloids [lj|, while the fission terms are empiri- 
cally defined from the typical behavior seen in the above 
simulations. The numbers 77,j change in time - we have 
{rij (i)}°^^, where rij (t) is the number of clusters of mass 
j at time t. 



(a) 





(b) 






FIG. 3: (a) The mean maximum cluster size M vs k for IBM 
simulations (A'' = 50). (b) p{m) as function of the cluster 
size m for rj = 0.34. Symbols show the average over eight 
IBM simulations for active particles with A^ = 50 and k = 1 
(circles) and k — 8 (crosses), errors bases give distributions 
of individual runs. The lines correspond to the mean field 
theory for k = 1 (solid) and k = 8 (dashed) (color online). 



This description neglects the geometry of clusters as 
well as spatial fluctuations. This allows us to consider 
a single rate constant for all possible collision processes 
between clusters of mass i and j , as well as a unique dis- 
integration constant for any cluster of mass i. In addition 
we make four crucial assumptions: i) The total number of 

particles in the system, N = J2i=i J''^j (^)i i^ conserved, 
ii) Only binary cluster collisions are considered. Colli- 
sions between any two clusters are allowed whenever the 
sum of the cluster masses is less or equal to N. iii) Clus- 
ters suffer spontaneous fission only by losing individual 
particles at the boundary one by one, i. e. a cluster can 
only decay by a process by which a j-cluster split into a 
single particle plus a (j — l)-cluster. This is motivated 
by observations in the above simulations, iv) All clusters 
move at constant speed, v « P/C\\ i which implies that 
rods in a cluster have high orientational order and in- 
teract only very weakly with their neighbors. Under all 
these assumptions the evolution of the n^'s is given by 
the following N equations: 
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where the dot denotes time derivative, Bj represents 
the fission rate of a cluster of mass j, defined by Bj = 
iv/R)^, and Aj^k is the collision rate between clusters 

of mass j and fc, defined by Aj^k = (va^i/A) ( ^/] + ^/k\ . 

(To is the scattering cross section of a single rod. R is 
the only free parameter and indicates the characteristic 
length a rod at the boundary of a cluster moves before 
it is leaving the cluster in a typical fission event. We as- 
sume R = aL taken into account that longer rods will 
stay attached to cluster for a longer time. 

Since ctq can be approximated hy ao ~ L + W — 
i/a ( y^ + -j= ] , the MFA depends only on the param- 
eters K, a. A, V and a. If one integrates Eqs. @ with 
parameters used in IBM simulations and an initial condi- 
0) = NSij, their solution yields steady state 
oo. From these values, we obtain a MFA 
for the weighted cluster size distribution p(m) — n^^ni/N 
for given values of the free parameters R resp. a. The 
best agreement between the MFA and the IBM simula- 
tions is found for a choice of a = 1.0±0.05 (see Fig. ^p). 
Hence, we will use i? = L in the following. To under- 
stand the relation between the parameters of the model 
and clustering effects, we can rescale Eqs. (0J by intro- 
ducing a new time variable: r = tv/ \/aK. The resulting 
equations *XU\ depend only on a dimensionless parameter 
P = [n + Vja/A. Note that w 7^ is scaled and does 
not affect the qualitative dynamics of the system. In the 
dimensionless model the parameter P stands for ratio be- 
tween fusion and fission processes and therefore triggers 
the transition from a unimodal to a bimodal cluster size 
distribution. We can easily establish a transition crite- 
rion, and by using a bisection method, we can accurately 
determine the critical transition parameter Pj.. Given the 
system area A, the rod area a and the number of rods 
iV, this method provides a way to calculate Kc'- 



K, - Pc{N)- - 1 
a 



(5) 
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At this point it is crucial that Pc depends on N , which 
is formally the number of equations in the MFA. By nu- 
merically solving the MFA equations for different par- 
ticle numbers iV up to TV = 1024, we find Pc{N) ex 
^-i.026±o.023_ rpj^jg j.ggyi|. indicates that in the MFA the 

critical parameter value Kc for the clustering transition 
does not depend on the number of particles. This does 
not imply that the weigthed distribution p(m) is inde- 
pendent of N\ in fact, we find that the probability for a 
rod to be in a large cluster increases with TV. We pro- 
ceed by assuming that Pc is inversely proportional with 
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FIG. 4: K-»7 phase diagram. The solid hne corresponds to the 
transition curve predicted by equation |S] The symbols indi- 
cate IBM simulations (A'' = 100). Crosses refer to unimodal 
pirn) and circles to bimodal p(r7i) of active particles. Stars 
refer to unimodal p(m) and hexagrams to bimodal •piri'i) of 
active-Brownian particles (color online). 



rods can achieve alignment at much lower densities than 
Brownian rods resp. particles in equilibrium systems. 
The clustering phenomenon is absent in simulations with 
isotropic self-propelled particles as well as with Brownian 
rods. Our model provides also an alternative explanation 
for collective behavior of rod-shaped objects - previous 
swarming models have achieved aggregation and cluster- 
ing by assuming attractive long-range interactions [J, la] ■ 
With respect to biology, our observation offers a simple 
physical explanation for the formation of clusters in many 
gliding rod-shaped bacteria, that often precedes the for- 
mation of biofilms and the appearance of more complex 
patterns. Acknowledgement: We acknowledge finan- 
cial support of Deutsche Forschungsgemeinschaft (DFG) 
through grant DE842/2 and fruitful discussions with L. 
Morelli, L. Brusch, J. Starruss and L. Sogaard- Andersen. 



A, we can express Kc as a simple function of the packing 
fraction: 



[1] 



Kc = CJT] ~ 1 



(6) 



where the constant is found to be C ~ 1.46. The k-t] 
phase diagram (Fig. Q shows a reasonable agreement of 
the transition line given by Eq. (jSJ and the IBM sim- 
ulation results. So, for the range of parameters used in 
the IBM, we retrieve in the MFA the unimodal shape of 
the weighted cluster size distribution for small values of 
K and 77, and the bimodal shape for large values of the 
two parameters. FiglJb gives a comparison of the cluster 
size distribution in the IBM and MFA. 

In summary, we have found non-equilibrium clustering 
for interacting self-propelled rod shaped-particles with 
sufhcient packing density 77 and aspect ratio k in sim- 
ulations. The rods interact via strong short range re- 
pulsive interactions that approximate excluded volume 
interactions. The onset of clustering has been defined 
by a transition from a unimodal to bimodal cluster size 
distribution. This transition is reproduced by a mean- 
field description of the cluster size distribution, which 
yielded a simple criterion, n = C/rj — 1, for the onset 
of clustering. This functional form with C ~ 1.46 pro- 
vides a good fit to the results of the simulations. The 
high density inside the cluster leads also to alignment of 
rods and coordinated motion of all particles in the clus- 
ter. The transition to clustering defined here is practi- 
cally independent of the system size resp. the number of 
particles. It is instructive to compare our result rewrit- 
ten in the form ktj + rj ~ 1.46 with the formula for the 
isotropic-nematic transition kij « 4.7 found in the two- 
dimensional version [l^ of Onsagers mean-field theory 
for Brownian rods loi. This shows that actively moving 
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